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Abstract 

We calculate the properties of an acoustic polaron in three dimensions in thermal 
equilibrium at a given low temperature using the path integral Monte Carlo 
method. The specialized numerical method used is described in full details, 
thus complementing our previous paper [R. Fantoni, Phys. Rev. B 86, 144304 
(2012)], and it appears to be the first time it has been used in this context. Our 
results are in favor of the presence of a phase transition from a localized state 
to an extended state for the electron as the phonon-electron coupling constant 
decreases. The phase transition manifests itself with a jump discontinuity in 
the potential energy as a function of the coupling constant and it affects the 
properties of the path of the electron in imaginary time: In the weak coupling 
regime the electron is in an extended state whereas in the strong coupling regime 
it is found in a self-trapped state. 

1. Introduction 

An electron in a ionic crystal polarizes the lattice in its neighborhood. An 
electron moving with its accompanying distortion of the lattice has sometimes 
been called a "polaron" [l], Q. Since 1933 Landau addresses the possibility 
whether an electron can be self-trapped (ST) in a deformable lattice Q [J 
This fundamental problem in solid state physics has been intensively studied for 
an optical polaron in an ionic crystal Q , Q, S, fl H, 111 | • Bogoliubov approached 



the polaron strong coupling limit with one of his canonical transformations. 
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Feynman used his path integral formalism and a variational principle to develop 



12j. Its extension 



an all coupling approximation for the polaron ground state 
to finite temperatures appeared first by Osaka [13I [l^ , and more recently by 
Castrigiano et al. [3, (3 12 1- Recently the polaron problem has gained new 
interest as it could play a role in explaining the properties of the high T c su- 
perconductors [lsj ]. The polaron problem has also been studied to describe an 
impurity in a Bose-Einstein ultracold quantum gas condensate of atoms 
In this context evidence for a transition between free and self-trapped optical 
polarons is found. For the solid state optical polaron no ST state has been found 
yet flflQ. 

The acoustic modes of lattice vibration are known to be responsible for the 

l|. Contrary to the optical mode which 



appearance of the ST state 
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interacts with the electron through Coulombic force and is dispersionless, the 
acoustic phonons have a linear dispersion coupled to the electron through a short 
range potential which is believed to play a crucial role in forming the ST state 
22 1 . Acoustic modes have also been widely studied Sumi and Toyozawa 



generalized the optical polaron model by including a coupling to the acoustic 
modes 2s| . Using Fcynman's variational approach, they found that the electron 
is ST with a very large effective mass as the acoustic coupling exceeds a critical 
value. Emin and Holstein also reached a similar conclusion within a scaling 
theory in which the Gaussian trial wave function is essentially identical to 
the harmonic trial action used in the Feynman's variational approach in the 
adiabatic limit 25 1. 



The ST state distinguishes itself from an extended state (ES) where the po- 
laron has lower mass and a bigger radius. A polaronic phase transition separates 
the two states with a breaking of translational symmetry in the ST one l| . The 
variational approach is unable to clearly assess the existence of the phase transi- 
tion [lj . In particular Gcrlach and Lowen Q concluded that no phase transition 
exists in a large class of polarons. The three dimensional acoustic polaron is 
not included in the class but Fisher et al. [2^] argued that its ground state is 
delocalized. 
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In a recent work 26j we employed for the first time a specialized path inte- 



gral Monte Carlo (PIMC) method 
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28( to the continuous, highly non-local, 



acoustic polaron problem at low temperature which is valid at all values of the 
coupling strength and solves the problem exactly (in a Monte Carlo sense) . The 
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method differs from previously employed methods 
hinges on the Levy construction and the multilevel Metropolis method with cor- 
related sampling. In such work the potential energy was calculated and it was 
shown that like the effective mass it usefully signals the transition between the 
ES and the ST state. Properties of ES and ST states were explicitly shown 
through the numerical simulation. 

Aim of the present paper is to give a detailed description of the PIMC 
method used in that calculation and some additional numerical results in order 
to complement the brief paper of Ref. [26]. 



The work is organized as follows: in section [5] we describe the acoustic po- 
laron model and Hamiltonian, in section [3] we describe the observables we are 
going to compute in the simulation, in section [4] we describe the PIMC numeri- 
cal scheme employed, in section[5]we describe the multilevel Metropolis method 
for sampling the path, in section[6]we describe the choice of the transition prob- 
ability and the level action, in section [7] we describe the correlated sampling. 
Section [5] is for the results, and section |H] is for final remarks. 



2. The model 

The acoustic polaron can be described by the following quasi-continuous 
model Q , 



2m 



k k 

Here x and p arc the electron coordinate and momentum operators respectively 
and a^, is the annihilation operator of the acoustic phonon with wave vector 
k. The first term in the Hamiltonian is the kinetic energy of the electron, 
the second term the energy of the phonons and the third term the coupling 
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energy between the electron and the phonons. The electron coordinate x is 
a continuous variable, while the phonons wave vector k is restricted by the 
Dcbyc cut-off k D . The acoustic phonons have a dispersion relation uj^ = uk 
(u being the sound velocity) and they interact with the electron of mass m 
through the interaction vertex Tk = huk (S / N) 1 / 2 [k / ko) 1 / 2 according to the 
deformation potential analysis of Rcf . {sfj . S is the coupling constant between 
the electron and the phonons and N the number of unit cells in the crystal with 
N/V = (47r/3)(fc /27r) 3 by Debyc approximation and V the crystal volume. 



Using the path integral representation (see Ref. [12| section 8.3), the phonon 
part in the Hamiltonian can be exactly integrated owing to its quadratic form 
in phonon coordinates, and one can write the partition function for a polaron 
in thermal equilibrium at an absolute temperature T (J3 = 1/ksT, with ks 
Boltzmann constant) as follows, 

r x=x(hp) 

ix=x(o) 



Z= dx e-T,s[x(t)Mt)A Vx{t) ; (2 ) 

J J-Jx=x(o) 





where the action S is given by 



37] 



= S,+U . (3) 

Here Sf is the free particle action, and IA the inter- action and we denoted with 

a dot a time derivative as usual. Using dimcnsionlcss units K = m = uk a = 

ks = V = 1 the action becomes, 

r 13 x 2 (t) r p f 3 
S= —^-dt + dt dsV eff (\x(t)-x(s)\,\t-s\) , (4) 
Jo * Jo Jo 

with the electron moving subject to an effective retarded potential, 

-r r S f i.f£q-(X(t)-X(s))-q\t-s\ /-n 

V *ff = l q<1 dqqe V ° (5) 



2I D 

35 ft 



\H \x(t)-x(s)\ L d qq 2 sin^a\x(t)-x(s)\ 



D -q\t-s\ 



(G) 



1 This is an approximation as e is neglected. The complete form is obtained by 

replacing e -"fcl*-s| by e - "*!' -3 /(I - e^ 1 "*) + e u * l t_s l e _ ' 3w '=/(l - e~' 9 "'=). But remember 
that f) is large. 
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where q = k/k , Id = f q <i = 47r/3, and we have introduced a non-adiabatic 
parameter 7 defined as the ratio of the average phonon energy, Huk to the 
electron band-width, (Kk ) 2 /2m. This parameter is of order of 10 -2 in typical 



ionic crystals with broad band so that the ST state is well-defined [23|. In 
our simulation we took 7 = 0.02. Note that the integral in ([6]) can be solved 
analytically and the resulting function tabulated. 

3. The observables 

In particular the internal energy E of the polaron is given by, 

idz 1 r , rr s ds /as\ 

where the internal energy tends to the ground state energy in the large (3 — > 00 
limit. 

Scaling the Euclidean time t = fit' and s = (3s' in Eq. ([¥]), deriving S with 
respect to /3, and undoing the scaling, we get, 

dS 1 f p x 2 , S [ p 7 

— = / — dt / dt dsx 

dP PJo 2 2I D J J 

dqqe^^-^-^^-qlt-sl) , (8) 

q<l P 

where the first term is the kinetic energy contribution to the internal energy, /C, 
and the last term is the potential energy contribution, V, 

3S f 3 f 1 fo(\/$9\*(t)-*( a )\) if 
p = -2Z dt ds dqq 3 L e -9\t-'\ x 



2/3 Jo Jo Jo J^q\x(t)-x(s 



So that, 



(2-g|i-*|) . (9) 



E = {K + V) . (10) 



An expression for K, not involving the polaron speed, can be obtained by 
taking the derivative with respect to (3 after having scaled both the time, as 
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before, and the coordinate x = \fj3x' '. Undoing the scaling in the end one gets, 

r/3 r 



K 



■J- f 'dt [ P ds [ dqqe ^i^m-x(s))- q \ t -s\ 

m D Jo Jo J g <i 



i\ -Q ■ (x{t) - x(s)) 
V 7 

3S f p , f , f 1 , „ 
-70 / * / ds / 

4 P Jo Jo Jo 

sin (^|g|aj(t) 



(11) 



cos ( \l —q\x(t) — x(s)\ 



$q\x(t) - x(a)\ 



,-9l*-*l 



(12) 



In the following we will explain how we calculated the potential energy P = (V). 
4. Discrete path integral expressions 

Generally we are interested in calculating the density matrix p = exp(—/3H) 
in the electron coordinate basis, namely, 



p(x a ,x b ;/3) = / / e~ b Vx{t) 
J-Jx=x„ 



(13) 



To calculate the path integral, we first choose a subset of all paths. To do 
this ,we divide the independent variable, Euclidean time, into steps of width 



r = p/M 



(14) 



This gives us a set of times, t% = kr spaced a distance r apart between and /? 
with fc = 0,l,2,...,M. 

At each time tk we select the special point Xk = x(tk), the k th time slice. 
We construct a path by connecting all points so selected by straight lines. It is 
possible to define a sum over all paths constructed in this manner by taking a 
multiple integral over all values of Xk for k = 1, 2, . . . , M — 1 where Xo = x a 
and xyi = Xb are the two fixed ends. The resulting equation is, 

• dxi dxM-i 



1 f f f -tdxi 

^,^)=M 3 j_J_ : ..^e- — ■ 

where the normalizing factor A = (2ttt) 3 ^ 2 . 



A 



(15) 



6 



The simplest discretized expression for the action can then be written as 
follows, 

s = j2 {Xk - 1 2T Xk) +T^Y.nu^) , (i6) 

k=l i=l j=l 

where V(ti,tj) = V e ff(\xi — Xj\, \i — j\) is a symmetric two variables function, 
V(s, t) = V(t, s). In our simulation we tabulated this function taking \xi— Xj\ = 
0,0.1, 0.2,..., 10 and \i - j\ = 0,1,..., M. 

In writing Eq. (fT6"|) we used the following approximate expressions, 

x k = Xk ~ Xk - 1 +0{t) , (17) 

T 

t k 



x\t)dt = ±iT + 0{T I ) , (18) 

fc-1 

/ I 3 V(s,t)dsdt = V{t i ,t j )T 2 + 0(t 3 ) . (19) 
Jti-i Jtj-i 

If we take V = in Eq. (fT6|) the M — 1 Gaussian integrals in (fT5|) can be done 
analytically. The result is the exact free particle density matrix, 

Pf (x a ,x b ;(3) = (2n(3)-^eM x °- x ^ . (20) 

Thus approximations (|17[) and (|18[) allow us to rewrite the polaron density 
matrix as follows, 

p(x a ,x b ;0) = J ■■ ■ J dx x ■ --dxM-i pf{x a ,xi;r) ■ ■ ■ p f (x M -i, x M ; t) x 

e-r'^-EjViuM . ( 2i) 

In the next section we will see that this expression offers a useful starting point 
for the construction of an algorithm for the sampling of the path: the Levy 
construction and the ana logy with classical polymer systems or the classical 



isomorphism described in 2jJ). 



The partition function is the trace of the density matrix, 

Z = J dxp(x,x;/3) . (22) 
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This restrict the path integral to an integral over closed paths only. In other 
words the paths we need to consider in calculating Z (and hence F) are closed 
by the periodic boundary condition, Xm = Xq = x. 

To calculate the internal energy we need then to perform the following M 
dimensional integral, 



oo poo 



E = — I I - I dxodxx-'-dxM^er 8 ^ + JC) 



(23) 

^J-ooJ-oo J-oo X M =X a 

To do this integral we use the Monte Carlo simulation technique described next. 



5. Sampling the path 

The total configuration space to be integrated over is made of elements 
s = [xq, Xi, . . . , Xm] where x^ are the path time slices subject to the periodic 
boundary condition xm = Xo- In the simulation we wish to sample these 
elements from the probability distribution, 

<s) = ^- , (24) 

where the partition function Z normalizes the function iz in this space. 

The idea is to find an efficient way to move the path in a random walk 
sampled by w, through configuration space. 

In order to be able to make the random walk diffuse fast thr oug h configura- 
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tion space, as r decreases, is necessary to use multislices moves 

In our simulation we chose to use the bisection method (a particular multi- 
level Monte Carlo sampling method j^)- That' s how an I levels move is con- 
structed. Clip out of the path m = 2 l subsequent time slices ajj, x^i, . . . , Xi+ m 
(choosing i randomly). In the first level we keep Xj and Xi +m fixed and, follow- 



ing Levy construction for a Brownian bridge 
at i + to/2 to, 



, we move the bisecting point 



Xi+m/2 = 7) \-V (25) 

where 77 is a normally distributed random vector with mean zero and standard 



deviation ^/tto/4. As shown in next section this kind of transition rule samples 
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the path using a transition probability distribution T oc exp(— Sf). Thus we 
will refer to it as free particle sampling. 

Having sampled x i+m / 2l we proceed to the second level bisecting the two 
new intervals (0, i + to/2) and (i + to/2, i + m) generating points x i+m /4 and 
:Ej_|_3 m /4 with the same algorithm. We continue recursively, doubling the number 
of sampled points at each level, stopping only when the time difference of the 
intervals is r. 

In this way we are able to partition the full configuration s into I levels, 
s = (s ,si, . . . ,si) where: s = [x , . • .,Xi,x i+m , . . .,x M -i], unchanged; Si = 
[fJi+m/a]) changed in level 1; s 2 = [x i+m/4> x i+3m/A }, changed in level 2; 
Si = [jEj+i, x i+2 , . . . , aj,- +m _i] changed in level I. 

To construct the random walk we use the multilevel Metropolis method 
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27j . Call (s[, . . . , s',) the new trial positions in the sense of a Metropolis 



rejection method, the unprimed ones are the corresponding old positions with 
s = s' . 

In order to decide if the sampling of the path should continue beyond level 
k, we need to construct the probability distribution itk for level k. This, usually 
called the level action, is a function of so, si . . . , Sk proportional to the reduced 
distribution function of s^ conditional on So,Si . . . , Sk-i- The optimal choice 
for the level action would thus be, 

7ifc(s , si . . . , sk) — J dsk+i ■ ■ ■ dsi 7r(s) . (26) 

This is only a guideline. Non optimal choices will lead to slower movement 
through configuration space. One needs to require only that feasible paths 
(closed ones) have non zero level action, and that the action at the last level be 
exact, 

iri(s , S!, . . . , si) = 7r(s) . (27) 

Given the level action 7Tfc(s) the optimal choice for the transition probability 
Tk(sk), for Sfe contingent on the levels already sampled, is given by, 

T ^) = ^% • (28) 
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One can show that Tt will be a normalized probability if and only if -n k is 
chosen as in (|26p . In general one need to require only that T k be a probability 
distribution non zero for feasible paths. In our simulation we used the free 
particle transition probability of the Levy construction as a starting point for a 
more efficient correlated sampling that will be described in a later section. 

Once the partitioning and the sampling rule are chosen, the sampling pro- 
ceeds past level k with probability, 

T k (s k )n k (s')TT k -i(s) 



A k (s ) = min 



1, 



(29) 



T k (s' k )ir k (s)ir k -i{s') _ 
That is we compare A k with a uniformly distributed random number in (0, 1), 
and if A k is larger, we go on to sample the next level. If A k is smaller, we make 
a new partitioning of the initial path, and start again from level 1. Here ttq 
needed in the first level can be set equal to 1, since it will cancel out of the 
ratio. 

This acceptance probability has been constructed so that it satisfies a form 
of "detailed balance" for each level k, 

^L Tk (s' k )A k (s') = ^l-T k (s k )A k (s) . (30) 

TTfc-l(s) 7Tfc_l(S'J 

The moves will always be accepted if the transition probabilities and level actions 
are set to their optimal values. 

The total transition probability for a trial move making it through all I levels 

is. 

i 

P(s^ s ') = l[T k ( s >)A k ( s i) . (31) 

k=l 

By multiplying Eq. (f3T)]) from k = 1 to k = I and using Eq. (f2~T|) . one can verify 
that the total move satisfy the detailed balance condition, 

tt(s)P(s s') = tt(s')P(s' s) . (32) 

Thus if there are no barriers or conserved quantities that restrict the walk to 
a subset of the full configuration space (i.e. assuming the random walk to be 
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ergodic) the algorithm will asymptotically converge to 7r, independent of the 
particular form chosen for the transition probabilities, T^, and the level actions, 

TTfc 



41 1 . We will call equilibration time the number of moves needed in the 
simulation to reach convergence. 

Whenever the last level is reached, one calculates the properties {JC and V) 
on the new path s', resets the initial path to the new path, and start a new 
move. We will call Monte Carlo step (MCS) any attempted move. 

6. Choice of and 

In our simulation we started moving the path with the Levy construction 
described in the preceding section. We will now show that this means that we 
are sampling an approximate T* with free particle sampling. 

For the free particle case (U — 0) one can find analytic expressions for the 
optimal level action 7r£ and the optimal transition rule T£. For examples for 
the first level, Eq. gives, 

n*(xi+ m /2) oc pf(xi,x i+m / 2 ;Tm/2)p f (x i+m /2,Xi +m ;Tm/2) (33) 
K e ^{Xi-x i+m/2 f e ^(x i+m/2 -x i+m ) 2 ( 34 ) 



x 



n/2 



X.+X, 



(35) 



This justify the Levy construction and shows that it exactly samples the free 
particle action (i.e. = 1 for all /c's). This also imply that for the interacting 
system we can introduce a level inter action, tt^ such that, 



TTfc = / dsk+i ...dsi 7r(s) , (36) 

with 



n{s) = — ■ (37) 



So that the acceptance probability will have the simplified expression, 

7rjb(s')#fc-i(s) 



Ak{s') = min 



1. 



TTk(s)TTk-l{s') 



(38) 
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For the k level inter action we chose the following expression, 



7Tfc oc exp 



[M/e k ] [M/e k ] 

«=1 j=l 



(39) 



where = m/2 k . In the last level £i = 1 and the level inter action 7r/ reduces 
to the exact inter action tt thus satisfying Eq. (|27|) . 

It' s important to notice that during the simulation we never need to calculate 
the complete level inter action since in the acceptance probabilities enter only 
ratios of level inter actions calculated on the old and on the new path. For 
example if for the move we clipped out the interval ti, . . . , tj+ ro with i + m < M 
[^|, we have, 



In 



T T 



-(r4) 2 < J2 v{ * u + m4r ' <i + n ^ r )+ 



m— n— 
Af 2 fe 



i-1 2 K 

5^ 51 + n£ k r) + Yl V ( mi kT, U + nt k r) ) 

m— i+m+1 n— 



m— 1 n— 



which is computationally much cheaper than ([39 



7. Correlated sampling 

When the path reaches equilibrium (i.e. P(s —> s') ss tt(s')) if we calculate, 



a(t / T ) 



aj(t) 



a;(i + i ) + sc(i-i ) 



(41) 



we see that these deviations are generally smaller than the free particle standard 
deviations used in the Levy construction (see Fig. [T]), 



2 When i + m > M there is a minor problem with the periodic boundary conditions and 
Eq. I|40j l will change. 
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Figure 1: Shows the deviations H4U for a simulation with S = 60 and S = 52.5, r = 0.025, 
I = 9. The free particle standard deviations 1 142 I I arc plotted for comparison. For S = 60 the 
path is localized while for S = 52.5 is unlocalized i.e. closer to the free particle path. 



(7/(4) - Vhr/2 



(42) 



As Fig. [T] shows, the discrepancy gets bigger as Ik increases. 

Wc thus corrected the sampling rule for the correct deviations. For example 
for the first level we used, 



T x {x i+m/2 )<xe 2 ^<-/ 2 > , 
where x = (xi + Xi+ m )/2. Since the level action is given by, 

( X i + rri/2- X ) 

iri(x i+m / 2 ) oc e 2 -/('™/ 2 > TTi{x l+m/2 ) , 



we can define a function, 



Pi oc e 



( X i + m /2- X ) 



[" 2 (™/2) <rJ(m/2)] 



and write the acceptance probability, 

P^s) 7ri(s')^o(s) 



Ai(s') = min 



1. 



Pi(s') 7ri(a)7f (a') 



(43) 



(44) 



(45) 



(46) 
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Which is a generalization of Eq. (|38p . 

We maintain the acceptance ratios in [0.15,0.65] by decreasing (or increas- 
ing) the number of levels in the multilevel algorithm as the acceptance ratios 
becomes too low (or too high). 

In the Appendix we report some remarks on the error analysis in our MC 
simulations. 

8. Numerical Results 

We simulated the acoustic polaron fixing the adiabatic coupling constant 
7 = 0.02 and the inverse temperature j3 = 15. Such temperature is found to 
be well suited to extract close to ground state properties of the polaron. The 
path was M time slices long and the time step was r = /3/M. For a given 
coupling constant S we computed the potential energy P extrapolating (with 
a linear \ square fit) to the continuum time limit, t — > 0, three points corre- 
sponding to time-steps choosen in the interval r E [1/100, 1/30]. An example 
of extrapolation is shown in Fig. [5] for the particular case j3 = 15, 7 = 0.02, and 
S = 60. 
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0.000 0.005 0.010 0.015 0.020 0.025 0.030 

x 

Figure 2: Shows the time step, r, extrapolation for the potential energy, P = (P). We run at 
/3 = 15,7 = 0.02, and S = 60. The extrapolated value to the continuum limit is in this case 
P = —16.1(5) which is in good agreement with the result of Ref. l33l . 

In Fig. [S] and Tab. [T] we show the results for the potential energy as a 
function of the coupling strength. With the coupling constant S = 52.5 we 
generated the equilibrium path which turns out to be unlocalized (see Fig. [4]). 
Changing the coupling constant to S = 60 and taking the unlocalized path as 
the initial path we sow the phase transition described in Fig. [3] the path after 
the phase transition is localized (see Fig. 01 . 
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Figure 3: At S = 60 the results for the potential energy V at each MC block (5 X 10 3 MCS) 
starting from an initial unlocalizcd path obtained by a previous simulation at S = 52.5. We 
can see that after about 30 blocks there is a transition from the ES state to the ST state. In 
the inset is shown the autocorrelation function, defined in Eq. (|A.8[I . for the potential energy, 
for the two states. The correlation time, in MC blocks, is shorter in the unlocalizcd phase 
than in the localized one. The computer time necessary to carry on a given number of Monte 
Carlo steps is longer for the unlocalizcd phase. 



ES 
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Figure 4: The top panel shows the polaron (closed) path x(t) as a function of Euclidean time 
t in units of r at equilibrium during the simulation. The middle panel shows the projection 
on the x — y plane of the path. The bottom panel shows the three-dimensional path. We see 
clearly how both path has moved from the initial path located on the origin but the path at 
S = 52.5 is much less localized than the one at S = 60. 
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Note that since S and r appear in the combination St 2 in U (and St in F) 
the same phase transition from an ES to a ST state will be observed increasing 
the temperature. With the same Hamiltonian we are able to describe two very 
different behaviors of the acoustic polaron as the temperature changes. 

In Fig. [5] we show the behavior of the potential energy as a function of 
the coupling strength. The numerical results suggests the existence of a phase 
transition between two different regimes which corresponds to the so called ES 
and ST states for the weak and strong coupling region respectively. We found 
that paths related to ES and ST are characteristically distinguishable. Two 
typical paths for the ES and ST regimes involved in Fig. [5] is illustrated in Fig. 
01 The path in ES state changes smoothly in a large time scale, whereas the path 
in ST state do so abruptly in a small time scale with a much smaller amplitude 
which is an indication that the polaron hardly moves. The local fluctuations 
in the results for the potential energy has an autocorrelation function (defined 
in Eq. (|A.8[1 ) which decay much more slowly in the ES state than in the ST 
state as shown in the inset of Fig. [3] Concerning the critical property of the 
transition between the ES and ST states our numerical results are in favor of 
the presence of a discontinuity in the potential energy. In the large j3 limit at 
/3 = 15 and fixing the adiabatic coupling constant to 7 = 0.02, the ST state 
appear at a value of the coupling constant between S = 52.5 and S = 55. With 
the increase of j3, the values for the potential energy P — (V) increase in the 
weak coupling regime but descrease in the strong coupling region. 

From second order perturbation theory (see Ref. 12j section 8.2) follows 
that the energy shift E(j, S) is given by — 3Sj[l/2 — 7 + 7 2 ln(l + I/7)] from 
which one extracts the potential energy shift by ta king P(j, S) = r ydE(j, S) j drf. 



231 follows that in the weak 



From the Feynman variational approach of Ref. 
regime the energy shift is —3S , 7[l/2— 7+7^(1 + 1/7)] and in the strong coupling 



regime —S + 3 1/5/57. 
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Figure 5: Shows the behavior of the potential energy P as a function of the coupling constant 
S. The points are the MC results (see Tab. [TJ, the dashed line is the second order perturbation 
theory result (perturbation) valid in the weak coupling regime and the dot-dashed line is the 
variational approach from Ref. [23l | (variational) in the weak and strong coupling regimes. 



9. Conclusions 

We used for the first time a specialized path integral Monte Carlo method 
to study the low temperature behavior of an acoustic polaron. At an inverse 
temperature /? = 15 (close to the ground state of the polaron) and at a non- 
adiabatic parameter 7 = 0.02 typical of ionic crystals we found numerical ev- 
idence for a phase transition between an extended state in the weak coupling 
regime and a self-trapped one in the strong coupling regime at a value of the 
phonons-electron coupling constant S = 54.3(7). The transition also appears 
looking at the potential energy as a function of the coupling constant where a 
jump discontinuity is observed. Comparison with the perturbation theory and 
the variational calculation of Ref. 23J is also presented. 



The specialized path integral Monte Carlo simulation method used as an 
unbiased way to study the properties of the acoustic polaron has been presented 
in full detail. It is based on the Levy construction and the multilevel Metropolis 
method with correlated sampling. Some remarks on the estimation of the errors 
in the Monte Carlo calculation are also given in the Appendix. This complement 
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Table 1: MC results for P as a function of S at /S = 15 and 7 = 0.02 displayed in Fig. [5] The 
runs where made of 5 X 10 5 MCS (with 5 X 10 4 MCS for the equilibration) for the ES states 
and 5 X 10 6 MCS (with 5 X 10 5 MCS for the equilibration) for the ST states. 



s 


P 


10 


-0.573(8) 


20 


-1.17(2) 


30 


-1.804(3) 


40 


-2.53(3) 


50 


-3.31(4) 


53.5 


-3.61(1) 


55 


-11.4(3) 


60 


-16.1(5) 


70 


-23.3(3) 


80 


-30.0(3) 



261 ] where fewer details on the Monte Carlo method had been 



our previous paper 
given. 

Our method differs from previously adopetd method 
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28 



29 



30 



31 



33, 



29J our path integral is in real space 
35| put the polaron on a lattice and 
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28| use PIMC single slice 



Unlike the method of Ref. 
rather than in Fourier space, Rcfs. [34] 
not on the continuum as we did, while Rcfs. 
move. The multilevel PIMC we used instead is a general sampling method 
which can efficiently make multislice moves. The efficiency £ (see Appendix) for 
the potential energy increases respect to the single slice sampling because the 
coarsest movements are sampled and rejected before the finer movements are 
even constructed. 

Although our results are of a numerical nature and we only probed the acous- 
tic polaron for one value of the non-adiabatic parameter 7 we think that the 
analysis support the existence of a localization phase transition as the phonons- 
clectron coupling constant S is increased at constant temperature or as the 
temperature is decreased at constant S. More so considering the fact that the 
introduction of a cut-off parameter have shown to work successfully in renor- 
malization treatments. 
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Appendix A. Estimating errors 

Since asymptotic convergence is guaranteed, the main issue is whether con- 
figuration space is explored thoroughly in a reasonable amount of computer 
time. Let us define a measure of the convergence rate and of the efficiency of a 
given random walk. This is needed to compare the efficiency of different transi- 
tion rules, to estimate how long the runs should be, and to calculate statistical 
errors. 

The rate of convergence is a function of the property being calculated. Let 
0(s) be a given property, and let its value at step k of the random walk be O k - 
Let the estimator for the mean and variance of a random walk with N MCS be, 

N-l 

0=<O k >=-^2°k , (A.l) 

k=0 

a 2 (G) =< (O k - Of > . (A.2) 
Then the estimator for the variance of the mean will be, 



° 2 (o) = <(^E<^E°) 2 > 

k k 



< (Ok - O) 2 > +2^T < (Ot - 0)(Oj 




2 



J2<(O i -0)(O j -0)> 




y ' i<j 

(A.7) 

The quantity ko is called the correlation time and can be calculated given the 
autocorrelation function for A k = Ok — O. The estimator for the autocorrelation 
function, Cfc, can be constructed observing that in the infinite random walk, 
< AiAj > has to be a function of \i — j\ only. Thus the estimator can be 
written, 

A A -, N-k 

_ < A A k > _ 1 ^ 

Ck - &HO) ~ (N - k)o*(0) ^ AnAl+k ■ (A - 8) 
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So that the estimator for the correlation time will be, 

2 N 

fco = l + -^(AT-fc)c fc . (A.9) 

k=l 

To determine the true statistical error in a random walk, one needs to estimate 
this correlation time. To do this, is very important that the total length of 
the random walk be much greater than ko . Otherwise the result and the error 
will be unreliable. Runs in which the number of steps N 3> ko are called well 
converged. 

The correlation time gives the average number of steps needed to decorrelate 
the property O. It will depend crucially on the transition rule and has a mini- 
mum value of 1 for the optimal rule (while cr(C) is independent of the sampling 
algorithm) . 

Normally the equilibration time will be at least as long as the equilibrium 
correlation time, but can be longer. Generally the equilibration time depends 
on the choice for the initial path. To lower this time is important to choose a 
physical initial path. Since the polaron system is isotropic, we chose the initial 
path with all time slices set to 0. 

The efficiency of a random walk procedure (for the property O) is defined 
as how quickly the error bars decrease as a function of the computer time, 
£e> = l/(J 2 (0)NT = l/a 2 (0)koT where T is the computer time per step. The 
efficiency depends not only on the algorithm but also on the computer and the 
implementation. 
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